Roton spectroscopy in a harmonically trapped dipolar Bose-Einstein condensate 
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We study a harmonically trapped Bose-Einstein condensate with dipole-dipole interactions in a regime where 
a roton spectrum emerges. We show that the roton spectrum is clearly revealed in the static and dynamic 
structure factors which can be measured using Bragg spectroscopy. We develop and validate a theory based on 
the local density approximation for the dynamic structure factor. 
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A variety of atomic species with large magnetic dipoles 
have been cooled to the quantum degenerate regime in which 
a Bose-Einstein condensate (BEC) forms |T|. The key feature 
of these gases is that the particles interact via a dipole-dipole 
interaction (DDI) that is long-ranged and anisotropic 0. A 
fascinating prediction is that when such a BEC is tightly con- 
fined in the direction that the dipoles are polarized, a roton- 
like spectrum can emerge (3). A number of theoretical pro- 
posals for detecting roton features have been made including 
sensitivity to external perturbations (H, depression in the crit- 
ical velocity 016]], and signatures in density fluctuations (7) 
(c.f. 0). However, to date there has been no experimental 
evidence for roton properties in dipolar BECs (9). 

In this paper we study a dipolar BEC confined in a quasi- 
two-dimensional (quasi-2D) harmonic trap. We vary contact 
and dipole interaction parameters over a wide range and char- 
acterize the emergence of a roton through the static and dy- 
namic structure factors. These quantities closely relate to the 
observable for Bragg spectroscopy ifTOlfTTIl . and thus are read- 
ily measured in experiments. We note that Bragg spectroscopy 
has emerged as a flexible tool for investigating ultra-cold 
gases and has been applied to resonant Bose [12 ] and Fermi 
lfl"3l gases, quasi- ID Bose gases lfT4ll . and vortices in BECs 
ifTSll . Recently the first application of Bragg spectroscopy to 
a dipolar BEC has been made fTEl in a nearly spherical trap, 
and used to demonstrate an anisotropic speed of sound. 

Our calculations for the structure factors are based on solv- 
ing the non-local Gross-Pitaevskii equation (GPE) for the con- 
densate and the Bogoliubov de-Gennes (BdG) equations for 
the quasi-particle excitations. Our calculations are fully three- 
dimensional (3D), i.e. we do not make the quasi-2D approx- 
imation in which an ansatz for the condensate shape in the 
tightly confined direction is assumed (this approximation has 
been shown to be surprisingly inaccurate in the regime of in- 
terest 1 17 ]). We finish by developing a local density approxi- 
mation (LDA) theory that provides a reasonably accurate de- 
scription of our full theory. We emphasize that the regime 
of our study is appropriate to current experiments with mag- 
netic dipoles. We refer the reader to Ref. fT8l for a discussion 
of the pure-2D regime that might be realized in future polar- 
molecule experiments. 

The DDI potential between a pair of dipolar atoms is 

TT ( s 3#dd 1 -3cos 2 <9 

^ dd(r) = ^r |,|3 » (1) 

where g^d = MoMm/3> Mm is the dipole strength, and is 



the angle between the dipole separation r and the polarization 
axis, which we take to be the z direction. The atoms also inter- 
act via a contact interaction of strength g = 4tt ah 2 /m with a 
the s-wave scattering length. We take the atoms to be confined 
in a cylindrically symmetric trap J7(x) = ^mu 2 (p 2 + \ 2 z 2 ) 
of aspect ratio X = uj z /uj p . Here we consider tight axial con- 
finement (uj z ^> uj p ) to produce a quasi-2D trap favorable for 
the emergence of rotons. 
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FIG. 1 . (color online) Stability phase diagram and related excitation 
properties of quasi-2D uniform dipolar BEC. White and light-grey 
regions indicate where the BEC is dynamically stable. In the light- 
grey region the spectrum has a roton minimum. In the dark-grey 
and black regions the system is dynamically unstable. This can arise 
from modes at zero momentum (phonon instability - black region) 
or finite momentum (roton instability - dark grey region) developing 
imaginary parts. Subplots A-D show cases of the spectrum ([3]), with 
the real (solid line) and imaginary (dashed line) parts shown. 

It is useful to review the properties of a homogeneous 
dipolar gas (i.e. with uj p = 0), which admits a simple ana- 
lytic treatment under the quasi-2D approximation (assuming 
a Gaussian mode structure along the tight direction). In this 
case the tight direction can be integrated out and the Fourier 
transform of the in-plane interaction potential is l3l 

VMM = g + SddF±(k p a z /y/2), (2) 

where F±(x) = 2 — 3^/irxe x2 erfc (x), g = g/V2ira z and 
Sdd = gdd/VZnciz are the 2D interaction parameters, with 
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a z = ^h/muj z the ^-confinement length scale. For a con- 
densate of areal density n the Bogoliubov dispersion relation 
for in-plane modes (i.e. no z excitation) is 



2e(k p )nV 2 r>(k p ), 



(3) 



where e(k p ) = h 2 k 2 p /2m. In Fig. [I] we show the generic fea- 
tures of the quasi-2D system as the dipolar and contact inter- 
action parameters are varied. Notably, the system can become 
unstable through a phonon or roton instability where k p — > 
(case B) or k p ~ l/a z (case D) modes, respectively, soften 
and develop imaginary eigenvalues. Within the stable region 
we have indicated a sub-region where the dispersion relation 
has a roton feature i.e. a local minimum at finite k p (case C). 

We now turn to our main concern, a fully 3D calculation of 
the trapped system. To do this we numerically solve the non- 
local GPE for the unit normalized condensate orbital ^o( x ) 
and chemical potential /i, and then diagonalize the BdG equa- 
tions for the quasi-particle excitations {uj(x),Vj(x)}, with 
respective energies ej lfT9l . Our numerical method, similar 
to Ref. l20lL utilizes the cylindrical symmetry by employ- 
ing a Fourier-Hankel representation to separate the GPE and 
BdG equations into a set of 2D problems specified by the an- 
gular quantum number m z = 0, ±1, ±2, . . . [e.g. Uj{x) = 
Uj (p, z)e %niz( f ) ]. We also use a cylindrically cutoff of the inter- 
action potential ([T]) to improve numerical convergence 1211 . 

The results we present focus on a trap with A = 40, al- 
though we find qualitatively similar behavior for A > 10. 
We have chosen A = 40 as being sufficiently tight for the 
roton to emerge at a reasonably large k value, yet is an as- 
pect ratio that is readily achievable in experiments (e.g. [22 ]). 
For convenience we introduce D = 3Ng^m/47rh 2 a p and 
C = Ngm/47rh 2 a p as dimensionless parameters for the dipo- 
lar and contact interactions, respectively, with a p = yJh/muo p 
and N the number of condensate atoms. 

In Fig. [2ja) we show the excitation spectrum of m z =0 
modes for C = as a function of D. For D > 400 we 
observe that high-lying quasi-particle modes begin to rapidly 
decrease in energy as D increases. We identify the soften- 
ing of these highly excited modes (i.e. modes with many ra- 
dial nodes) as the manifestation of the roton spectrum in the 
trapped gas. This interpretation is supported by other stud- 
ies of the trapped system in which the quasi-particle spec- 
trum was approximately mapped onto a dispersion relation 
l23l (also see k defined below). At D « 730 the first of these 
roton (quasi-particle) modes hits zero energy and develops an 
imaginary part, signaling the onset of a dynamical instability. 
We note that modes with \m z \ > exhibit similar trends to 
the m z = results shown in Fig.[2ja). 

The condensate orbital and the density perturbations associ- 
ated with two quasi-particle modes are indicated in Fig. [2jb)- 
(d). Notably the roton mode [Fig. [2jd)] is localized near the 
center of the condensate and has a short wavelength. Follow- 
ing l23l we assign a wavevector to the modes according to 

k = yj (k 2 ), and find that for this mode k = 6A5/a p , similar 

to the inverse z confinement length l/a z w 6.3/ a p . 
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FIG. 2. (color online) (a) Spectrum of m z = quasi-particle modes 
demonstrating the roton- softening, (b)-(d) Contour plots of modes at 
D — 700. (b) Condensate density, (c) and (d) give the density per- 
turbation [Spj =ipo(uj+Vj)] for the first [(c)] and third (roton) [(d)] 
m z = quasi-particles [these modes are indicated by circles in (a)]. 
Solid (dotted) contours indicate positive (negative) density perturba- 
tions, (e) S(k) for D values corresponding to the vertical dashed 
lines in (a). Inset: A comparison of full calculations for S(k) (sym- 
bols) to the LDA results [Eq. f6}] (lines) for D = 200, 600, 725 
(lowest to higher curves). Results for A = 40 and C = 0. 



The dynamic structure factor for a T = BEC is ifTTl . 



S(k,o;) = 



dxK(x) +^*(x)K k - x ^o(x) S(uj-uJj) 



(4) 

where ujj = ej/h. It is worth noting that Bragg spec- 
troscopy measures the imaginary part of the response func- 
tion Im[xkM] = -|[S(k,o;) - 5(-k, -uj)} EH, and to 
leading order this is only sensitive to the zero-temperature dy- 
namic structure factor. Thus our results should be applicable 
to regimes with a discernible non-condensate fraction. Cor- 
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rections beyond leading order will require a finite temperature 
extension of the theory (e.g. see l25l ). 

Integrating S(k,uj) over frequency yields the static struc- 
ture factor S(k) = J dujS(k,uj), which also relates to the 
Fourier transform of the pair correlation function [11]. For 
the uniform system S'(k) directly gives the dispersion relation 
through the Bijl-Feynman formula S(k) = e(k)/e£(k). 

Here we restrict our attention to evaluating the structure fac- 
tors for radial wavevectors k = k k p , since the roton modes 
exhibit non-trivial structure in-plane [see Fig.[2jd)], and from 
hereon will denote these with scalar arguments, i.e. S(k,u), 
S(k). For a given value of k the numerical evaluation of 
S(k,uj) requires including modes up to a maximum energy 
e max with e max > h 2 k 2 /2m. In practice we check that suffi- 
ciently many modes are included by ensuring that the /-sum 
rule, J °° duj(jjS(k, uj) = Hk 2 /2m, is satisfied. For the k val- 
ues we consider here (k < 20/a p ) we typically use > 10 4 
modes in our calculations. 

We present results for S(k) in Fig.|2je) for various values 
of the dipole interaction strength. The suppression of S(k) as 
k —> reveals the low-energy phonon spectrum of the sys- 
tem (see ifTOl ). A significant peak in S(k) at /c pea k ~ 6.5/a p 
forms for interaction values of D > 400, which corresponds 
to where the high energy modes begin to rapidly descend in 
the spectrum [see Fig.[2ja)]. Appealing to the Bijl-Feynman 
formula we identify a significant peak in the static structure 
factor with the appearance of a roton feature in the excitation 
spectrum. This identification is useful because it corresponds 
to a practical experimental observable and does not depend 
upon any ad hoc scheme for assigning a dispersion relation to 
the excitations of the trapped system. 
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FIG. 3. (color online) Characterization of the roton properties of a 
A = 40 trapped dipolar BEC using S(k). (a) Peak value of S(k). 
(b) Peak wavevector /c pea k (none shown when max[5'(/c)] < 1.05). 
Circles mark the parameters of states analyzed in Fig. |2je). Grey 
boundary line indicates when the system is unstable due to excited 
modes softening, with black segments indicating that the BEC is in a 
bi-concave state at the boundary (see 0). 

In Fig. [3] we characterize the behavior of S(k) over a broad 



range of contact and dipole parameters where the BEC is dy- 
namically stable. We show where a peak in S(k) emerges 
and characterize its height [Fig.[3ja)] and wavevector (fc pea k) 
[Fig. [3jb)]. Our results show that the roton character of the 
spectrum is generally enhanced [i.e. height of peak in S(k) 
increases] at fixed dipole strength by decreasing (i.e. making 
more negative) the contact interaction strength, although the 
value of fcpeak tends to decrease as this happens. 

In Fig. [3] we also indicate the boundary upon which the 
BEC becomes dynamically unstable. Because the cloud is 
tightly confined in the z direction the repulsive character of 
the DDI (due to side-by- side dipoles) dominates and hence 
the DDI partially stabilizes the BEC against collapse from 
a negative value of the contact interaction. Similar observa- 
tions, for gases with smaller trap aspect ratio, were presented 
in Ref. l2Tll . There are regions near the boundary where the 
condensate develops a bi-concave density profile with a local 
minimum in the BEC density at trap centre [3 ] [26 ]. We do 
notice any signature of the bi-concave BEC in S(k). 
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FIG. 4. (color online) (a)-(f) S(k, uj) for indicated values of D. A 
Gaussian of width Acj = 0.5 uj p is used to smooth the ^-functions 
in S(k, uj). The white dashed line shows the free particle dispersion 
e(k)/h. The solid black and grey lines are the mean response Cu(k) 
obtained from GPE based and quasi-2D approximation LDA calcula- 
tions, respectively [see text and below]. White ellipse in (f) identifies 
a roton response feature. LDA calculations of S(k, uj) for the param- 
eters in (f) are shown in (fl) and (f2). (fl) GPE based LDA approach 
of Eq. (6). (f2) quasi-2D approximation LDA approach (see text). 
The mean response u(k) = f duo ujS(k,uj)/S(k), is shown for each 
result. Other parameters: C = and A = 40. 

Because the frequency dependence of the system response 
is most directly measured in Bragg spectroscopy experiments 
it is also worth discussing the behavior of the dynamic struc- 
ture factor. In Figs. |5Ja) -(f) we show S(k,uj) for the same 
cases considered in Fig. [2je). In the vicinity of the ro- 
ton wavevector [i.e. k ~ 6.5/a p ] the frequency response is 
quite broad and dips down sharply towards zero frequency for 
D > 700. The discernible response feature indicated with 
an ellipse in Fig. |4jf) is due to the roton mode identified in 
Fig. |2jd) (but also has contributions from similar modes with 

Kl >o). 
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For the case of BECs with contact interactions a success- 
ful analytic approximation for S (k, uj ) has been developed us- 
ing the Thomas-Fermi approximation for the condensate and 
treating the excitations within LDA ifTTIl . We have found that 
the main issue in extending this type of analysis to the tightly 
confined dipolar BEC arises from the sensitive dependence 
of the in-plane DDI potential [c.f. Eq. ([2])] upon the shape of 
the condensate in the z direction. For this reason we use the 
GPE solution itself as the basis for calculating S(k,u) using 
an LDA treatment of the excitations, thus avoiding the need to 
numerically solve for the BdG equations. We note that gener- 
alized z mode treatments, e.g. Ref. IfTTIl . could also be used. 
We define a locally varying in-plane interaction potential 



dk z 



27rn 2 (p) 



+ Udd(k p ,k z ) [h(k z ,p)} 2 , (5) 



where n{p) = j dz \ipo(p, z)\ 2 is the areal density, h(k Z: p) 
is the z-Fourier transform of the condensate density, and 
Udd(k p ,k z ) = gd(3kl/k 2 — 1) is the Fourier transform of 

£/dd( r )- The p dependence of V^ D (/c p ,p) accounts for the 
changing z profile of the condensate as p varies. Note that 
in the limit of vanishing interactions, where the z shape of the 
condensate is a Gaussian (independent of p), V^ D (/c p ,p) re- 
duces to the analytic result in Eq. We construct S(k,uj) 
treating the in-plane excitations with the LDA, i.e. summing 
over the parts of the BEC at various densities IfTTIl 

Slda(^)= / dp 2np n ^ k p] 5 (uj-e B (k p , p)/h) , (6) 
J £B{k p ,p) 



where e B (k p ,p) = y e(k p ) 2 + 2e(k p )n(p)V^ D (k p , p). In 
Figs. [4] (fl) and (f2) we compare our GPE based LDA ([6]) 
against LDA calculation using the quasi-2D approximation 
[this only differs by the replacement V^j^fcp, p) ~~ ^ ^2d(^p) 
m Eq. {6}]. This comparison reveals the sensitivity of S(k, uj) 
to the z shape of the condensate. In the inset to Fig. [2je) we 
compare the full numerical calculations of S(k) against the 
GPE based LDA, and find good agreements until D > 700 
where the roton modes approach zero energy. 

Finally, we relate the dimensionless parameters of our cal- 
culations to current experimental systems. A value of D ~ 
700 for uj p = 2tt x 40 s _1 would require a condensate with 
N = {64.3, 4.1, 8.1} x 10 4 atoms for { 52 Cr, 164 Dy, 168 Er}, 
respectively. The value of D can be adjusted by changing the 
radial confinement, atom number or dipolar strength [ 27 ]. Our 
results demonstrate [see Fig.[3ja)] that instead, for fixed dipole 
strength, the roton spectrum can be accessed by making the s- 
wave scattering length negative (e.g. see (28)). 

In conclusion we have explored the excitation properties of 
a quasi-2D dipolar BEC in terms of the dynamic and static 
structure factors. Our results show that clear and direct signa- 
tures for the roton spectrum will emerge in the structure fac- 
tors and should be readily observable with Bragg spectroscopy 
in current experiments. We have constructed an approximate 
LDA theory for S(k, uj) which we have validated against the 
full theory. Future work will consider the extension to a set of 
quasi-2D traps realized with an optical lattice lfT71[29ll . 
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